/*! 
 ***************************************************************************
 *
 *  \file test2D.cpp
 *  \brief Particle Tracking classes: examples of use (2D)
 *
 *  \author  Luis M. de la Cruz Salas
 *  \date    Last Revision: Wed Set 13 2006
 *  \version 1.0
 *
 ***************************************************************************
 *
 *  \par Description
 *  This file presents some examples for using ParticleTracking<>, 
 *  Euler<>, RungeKutta4<> and Bilinear classes. These classes are intended 
 *  to handle some algorithms to track mass less particles inside flows 
 *  defined in rectangular and structured meshes. We work in two dimensions 
 *  in this example. Modify the PARTRACK and TOONA variables in the Makefile 
 *  to the right places and then type the next commands:
 *
    \verbatim
    % make test2D   (to compile)
    % ./test2D      (to run) \endverbatim
 *  \par
 *
 ***************************************************************************
 */

#include <iostream>
#include "Meshes/StructuredMesh.hpp"
#include "Meshes/Uniform.hpp"
#include "Utils/inout.hpp" 

#include "ParTrack/ParticleTracking.hpp"
#include "ParTrack/Linear.hpp"
#include "ParTrack/Euler.hpp"
#include "ParTrack/RungeKutta4.hpp"

using namespace std;
using namespace Tuna;

typedef TunaArray<double, 2>::tiny floatTinyArray_t;
typedef TunaArray<double, 2>::huge ScalarField2D;

int main () 
{
    double dt = 0.0001;
    floatTinyArray_t positionE, positionR;
    positionE(0) = 0.3;
    positionE(1) = 0.2;
    positionR = positionE;

    StructuredMesh< Uniform<double, 2> > mesh(1, 256, 1, 256);
    floatTinyArray_t deltas = mesh.getDeltas();

    mesh.print();

    cout << "\n ----- Tracking: begins -----\n";
    cout << " Time step   : " << dt << endl;
    cout << " Initial pos : " << positionE << endl;
    cout << " Mesh size   : " << deltas << endl;    

    cout << "\n ----- Tracking: Euler Method -----\n";

    ScalarField2D u(mesh.getExtentNodes());       
    ScalarField2D v(mesh.getExtentNodes());
    InOut::readFromFile_1(u, v, 151, "./Velocity/velc.");
    
    ParticleTracking< Euler< Linear <double, 2> > > particleE;
    positionE = particleE.track(positionE, deltas, u, v, dt);
    cout << " New pos : " << positionE << endl;

    cout << "\n ----- Tracking: Runge-Kutta 4to order Method -----\n";
//
// Read the velocity field at time t + dt
//
    ScalarField2D un(mesh.getExtentNodes());       
    ScalarField2D vn(mesh.getExtentNodes());
    InOut::readFromFile_1(un, vn, 152, "./Velocity/velc.");

    ParticleTracking< RungeKutta4< Linear<double, 2> > > particleR;
    positionR = particleR.track(positionR, deltas, u, v, un, vn, dt);
    cout << " New pos : " << positionR << endl;

    cout << "\n ----- Tracking: ends -----\n";

    return 0;
}



